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Characterization of Extinction and Reignition of 
Hydrogen- Air Diffusion Flame 


Abstract 

Work performed during the first six months of the project duration for i 
the NASA Grant (NAG-1-861) is reported here. An analytical and computa- I 

tional study of opposed jet diffusion flame for the purpose of understand- ' 

| 

ing the effects of contaminants in the reactants and thermal diffusion of I 
light species on extinction and reignition of diffusion flames is in prog- I 
ress. The methodologies that have been attempted so far are described. I 

Introduction 

Extinction and reignition of diffusion flame are of major concern to the 
designers of aerospace propulsion systems. Engine reliability can be improved 
by knowing under what conditions extinction and reignition takes place. Velo- 
city gradient, turbulence, presence of contaminants, and characteristics of 
fuel are some of the factors that affect the above mentioned phenomena. Sev- 
eral attempts have been made in the past to study these phenomena by exper- 
imental as well as analytical means . One commonly used experimental setup is 
called the opposed jet burner (Figure 1). For purposes of discussion, an axi- 
symmetric configuration will be used here. Fuel and oxident gases, often 
diluted with nitrogen, flow out of axially aligned tubes of equal diameter; 
the two gaseous streams flow towards each other, and finally meet, creating a 
stagnation point flow similar to the classical stagnation point flow encoun- 
tered in boundary- layer theory. Under favorable conditions, a flat diffusion 
flame is formed near the stagnation plane. The flame may be made to blow off 
or reignite in the stagnation region by changing the parameters such as flow 
velocity, fuel/oxident dilution level, tube diameter, and tube separation dis- 


tance . 



Extensive experiments have been conducted at NASA Langley Research Center 
on the opposed jet burner. Effects of the presence of contaminants, fuel 
dilution level, velocity gradient, tube size and tube separation have been 
extensively examined. A hysteresis phenomenon in the context of extinction 
and reignition, previously alluded to in the literature, has been observed and 
systematically documented for the first time. 


Analysis and Numerical Solution 

The governing equations of continuity, momentum, energy, and species con- 
tinuity are transformed into a set of coupled ordinary differential equations 
using the boundary -layer approximations and the Howarth-Dorodnytzyn transfor- 
mation. The resulting set of equations are given below. 
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In the above expressions C — pp/PcoPco* Pr “ Prandtl number, Le - Lewis number, 
Y - mass fraction, <!>' * ' - species production rate (mass/unit time/unit vol- 
ume) . Subscript i denotes the ith species. 

The problem is thus reduced to a two -point boundary -value problem. The 
boundary conditions are as follows. 
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Several efficient routines are available for the solution of multi-point 
boundary -value problems. The routine COLSYS, obtained from the mathematical 
library at LaRC f has been found to work well for the present problem and, 
therefore, it will be used for the present study. 

Chemical Kinetics 

The nonequilibrium nature of the flame requires chemical kinetics infor- 
mation for the calculation of chemical species production rates. A library of 
routines known as CHEMKIN, capable of generating chemical kinetics informa- 
tion, and thermodynamic and transport properties, has been developed at Sandia 
National Labs. These routines have been developed in modular form and are 
totally problem- independent . They have already been tested and are found to 
work well for the present problem. The principal investigator has already 
incorporated the CHEMKIN routines into the computer program developed for the 
present problem. 


Work Completed 

Numerical algorithms for the governing equations described earlier have 
been developed and solution attempted on the NASA Langley computer system 
as well as the computer system at the University of Missouri-Rolla. The 
following are brief descriptions of the algorithms. 

1. Collocation Method: This is a routine called COLSYS [12] available 

from the Langley Math Library. The routine can be used to solve boundary- 
value problems governed by a system of ordinary differential equations. 
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2 . Fully Implicit Finite Difference Method: Recognizing the stiff nature 

of the differential equations because of the presence of the chemical 
source term, a fully implicit algorithm was written to avoid the limitation 
on step size for stable solution. 

3. Patankar and Spalding Method: This method uses the conservative form 

of the equations. The source term is treated according to the rules for- 
mulated by Patankar and Spalding [13]. 

Work in Progress 

Rogers and Chinitz [14] global chemistry model has been used in many 
situations for the solution of the hydrgen-air problem. According to this 
model, two global reactions are used to represent the combustion of 
hydrogen in air. Algebraic or differential equations based on the size of 
the source term, are used in the numerical algorithm. 
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Figure 2. Program Flow Chart for Analyzing Opposed 
Jet Diffusion Flame. 
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Appendix 


Thermal Diffusion Ratios 

The following Information is extracted from Reference 11. We consider 
that a thermal diffusion ratio 0 k can be defined such that the thermal diffu- 
sion velocity W k ^ is given by 

D k*k 1 3T 
X k T dXi 


Wki - 


(4) 


in which is a spatial coordinate. The mole fractions are given by X^, and 
the are "mixture diffusion coefficients" that are defined in a following 
section on average properties. We only consider thermal diffusion in the 
trace, light component limit (specifically, species k having molecular weight 
less than 5). The thermal diffusion ratio is given by 


in which 
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In the computer code, we store fits to j , B^j , and C^j , instead of evaluat- 
ing them from the preceeding three equations. We fit the data presented in 
tables by Monchick and Mason [7]. 
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The Fitting Procedure 


To expedite the evaluation of transport properties in a computer code, 
such as a flame code, we fit the temperature -dependent parts of the property 
expressions. Then, rather than evaluate the complex expressions for the 
properties, we need only evaluate the comparatively simple fits. The fitting 
procedure must be carried out for the particular system of gases that are 
present in a given problem. Therefore, the fitting can not be done "once and 
for all," but must be done once at the beginning of each new problem. 

We use a polynomial fit of the logarithm of the property versus the 
logarithm of the temperature. For example, for the viscosity, we have 

N 

ln, 7k “ I a n,k( lnT ) n ’ 1 < < 10 ) 

n=l 

and similarly for the thermal conductivity we have 

N 

lnA k - l bnkClnT)"- 1 . (11) 

n=l 

The fits are done for each pair of binary diffusion coefficients in the sys- 
tem. 

N 

lnDjk - X ^n,jk^ n ^) n (12) 

n~l 

The viscosities and thermal conductivities are independent of pressure, 
but the diffusion coefficients depend inversely on pressure. We compute all 
the coefficients at one atmosphere pressure. When later evaluating a diffu- 
sion coefficient, we simply divide the diffusion coefficient as evaluated from 
the fit by the actual pressure in atmospheres. 

The thermal diffusion ratios are weakly dependent on temperature. In 
fitting these we use polynomials in temperature, rather than logarithm tem- 
perature. Evaluating these polynomials is more computationally efficient. 
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Average Mixture Properties 


As the fits are done for the pure species viscosities, conductivities, 
and for the binary diffusion coefficients, we need a means to determine the 
properties of a gas mixture. We provide subroutines for the evaluation of 
mixture viscosities and conductivities via two different approaches. 

The first approach, and the simplest, is the combination averaging formu- 
las recommended by Mathur, et al. ® for the evaluation of conductivity. 


A 




(13) 


in which are the species mole fractions. The analogous formula for viscosity 
is, 


*7 
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The semi -empirical formulas due to Wilke and modified by others can also 
be used to evaluate mixture properties, but they are significantly more expen- 
sive computationally, and apparently more accurate. The Wilke formulas for 
mixture conductivity are given by 
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As the appears in the "inner loop", the computational expense is obvi- 

ous . 


There is a noticeable difference between the properties as evaluated by 
the two methods. At high temperature we find that the conductivities evalu- 
ated by the Wilke formulas are consistently higher (on the order of 10 percent 
at 2000K for hydrogen- air mixtures) than those evaluated by the combination 
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formulas. There is a much larger discrepancy between the two formulas for 
viscosities (sometimes as much as a factor of two). At least for viscosities, 
the Wilke formulas should be the more accurate. 

We assume that a diffusion coefficient for the diffusion of species k into 
the gas mixture can be defined such that the species diffusion velocity is 
related to the mole fraction gradient in the following way: 


A ki “ -°k 


l_«k 


(17) 


We determine the diffusion coefficient of the species k in the mixture by the 
following formula: 


i - Vk 
£j A ^ 


(18) 


in which are the species mass fractions, and the mole fractions. A poten- 
tial problem with this expression is that it is not mathematically well def- 
ined in the limit of the mixture becoming a pure species. Even though diffu- 
sion itself has no real meaning in the case of a pure species, we need to be 
sure that in a programming implementation the diffusion coefficients behave 
reasonably and that the code does not "blow up" when the pure species condi- 
tion is reached. We circumvent these problems by evaluating the diffusion 
coefficients in the following way. First, write 


°k 


1 - Y k 


K 
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and, as the mole fractions and mass fractions are related by 

x k w k 
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we can write 


x j w j 
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In this form the roundoff is accumulated in roughly the same way in both the 
numerator and denominator, and thus the quotient is well behaved as the pure 
species limit is approached. However, when we have exactly a pure species, 
the formula is still undefined. To overcome this difficulty we add a small 
perturbation to the mass fractions when computing the mole fractions. In 
other words, for the purposes of computing mixture diffusion coefficients, we 
simply do not allow a pure species situation to occur; we always maintain a 
residual amount of each species. Specifically, we assume in the above formu- 
las that 

WY k 

X k + 6 . (22) 

k % 

A value of 10 "^ for S works well. It is small enough that it is numerically 
insignificant compared to any mass fraction of interest, yet there is no 
trouble representing it on any computer. 

Species Conservation 

Some care needs to be taken in using the mixture diffusion coefficients as 
described here. The formulas are approximations, and they are not constrained 
to require that the net species diffusion flux is zero, i.e., the condition 

K 

l V ki Y k - 0 (23) 

k«l 

is not satisfied. Therefore, one must expect that applying these mtixture 
diffusion relationships in the solution of a system of species conservation 
equations should lead to some nonconservation, i.e., the resultant mass frac- 
tions will not sum to one. Therefore, one of a number of corrective actions 
must be invoked to ensure mass conservation. 

Unfortunately, resolution of the conservation problem requires knowledge 
of species flux, and hence details of the specific problem and discretization 
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method. Therefore, it is not reasonable in the general setting to attempt to 
enforce conservation. 

One attractive method is to define a "conservation diffusion velocity" as 
Coffee and Heimerl recommend. In this approach we assume that the dif- 
fusion velocity vector is given as 

V ki “ *ki + W kl + v cl , (25) 

in which i/ k ^ is the ordinary diffusion velocity due to mole fraction gra- 
dients, W k £ is the diffusion velocity due to temperature gradients, and V c ^ is 
a constant correction factor (independent of species, but spatially varying) 
introduced to satisfy Equation 24. The correction velocity is defined by 

K 

V ci “ ‘ I Y k (., ki + W ki ). (26) 

k=l 

This approach is the one followed by Miller, et al. in their flame models. 

An alternative approach is attractive in problems having one species which 
is always present in excess. Here, we do not solve a conservation equation 
for the one excess species. The excess component's mass fraction is then 
simply computed by subtracting the sum of the remaining mass fractions from 
one. A similar approach involves determining locally at each computational 
cell which species is in excess. The diffusion velocity for that species is 
computed to require satisfaction of Equation 23. 
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